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Abstract A new general relativistic magnetohydrody- 
namics (GRMHD) code "RAISHIN" used to simulate 
jet generation by rotating and non-rotating black holes 
with a geometrically thin Keplarian accretion disk finds 
that the jet develops a spine-sheath structure in the 
rotating black hole case. Spine-sheath structure and 
strong magnetic fields significantly modify the Kelvin- 
Hclmholtz (KH) velocity shear driven instability. The 
RAISHIN code has been used in its relativistic magne- 
tohydrodynamic (RMHD) configuration to study the 
effects of strong magnetic fields and weakly relativistic 
sheath motion, c/2, on the KH instability associated 
with a relativistic, 7 ~ 2.5, jet spine-sheath interac- 
tion. In the simulations sound speeds up to ~ c/VS 
and Alfven wave speeds up to 0.56 c are considered. 
Numerical simulation results are compared to theoreti- 
cal predictions from a new normal mode analysis of the 
RMHD equations. Increased stability of a weakly mag- 
netized system resulting from c/2 sheath speeds and 
stabilization of a strongly magnetized system resulting 
from c/2 sheath speeds is found. 

Keyvirords galaxies: jets — gamma rays: bursts — 
ISM: jets and outflows — methods: analytical — MHD 
— relativity — instabilities 

1 Introduction 

Relativistic jets are associated with active galactic nu- 
clei and quasars (AGN), with black hole binary sys- 
tems (microquasars), and are thought responsible for 
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the gamma-ray bursts (GRBs). The observed proper 
motions in AGN and microquasar jets imply speeds 
from ~ 0.9 c (e.g., Mirabel & Rodriquez 1999) up to 
- 0.999 c (e.g., the 3C345 jet Zensus et al. 1995; Stef- 
fen et al. 1995), and the inferred speeds for GRBs are 
~ 0.99999 c (e.g., Piran 2005). 

Jets at the larger scales may be kinetically domi- 
nated and contain relatively weak magnetic fields, but 
stronger magnetic fields exist closer to the acceleration 
and coUimation region. Here GRMHD simulations of 
jet formation (e.g., Koide et al. 2000; Nishikawa et 
al. 2005; De Villiers et al. 2003, 2005; Hawley & Kro- 
hk 2006; McKinney & Gammie 2004; McKinney 2006; 
Mizuno et al. 2006) and earlier theoretical work (e.g., 
Lovelace 1976; Blandford 1976; Blandford & Znajek 
1977; Blandford & Payne 1982) invoke strong magnetic 
fields. Additionally, the GRMHD simulations suggest 
that jets driven by magnetic fields threading the er- 
gosphere can reside within a broader sheath outflow 
driven by the magnetic fields anchored in the accre- 
tion disk (e.g., McKinney 2006; Hawley & Krolik 2006; 
Mizuno et al. 2006), or less coUimated accretion disk 
wind (e.g., Nishikawa et al. 2005). 

Recent observations of QSO winds with speeds, 
0.1 — 0.4c, also indicate that a jet could reside in a high 
speed sheath (Chartas et al. 2002, 2003; Pounds et al. 
2003a, 2003b; Reeves et al. 2003). Circumstantial evi- 
dence such as the requirement for large Lorentz factors 
suggested by the TeV BL Lacs when contrasted with 
much slower observed motions has been used to suggest 
the presence of a spine-sheath morphology (Ghisellini 
et al. 2005), and Siemignowska et al. (2007) have pro- 
posed a spine-sheath model for the PKS 1127-145 jet. 
Spine-sheath structure has also been proposed based 
on theoretical arguments (e.g., Sol et al. 1989; Henri & 
Pelletier 1991; Laing 1996; Meier 2003) and has been 
investigated in the context of GRB jets (e.g., Rossi et 
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al. 2002; Lazzatti & Begelman 2005; Zhang et al. 2003, 
2004; Morsony et al. 2006). 

In §2 we illustrate the spine-sheath configuration 
found by our GRMHD jet generation simulations. Pre- 
vious relativistic fluid dynamical (RHD) simulation and 
theoretical work has shown the importance of spine- 
sheath structure to KH instability (Hardee & Hughes 
2003). In §3 we report on numerical results that extend 
this previous investigation numerically and in §4 theo- 
retically to the strongly magnetized RMHD regime. 

2 GRMHD Jet Spine-Sheath Generation 

In order to study the formation of relativistic jets 
from a geometrically thin Keplerian disk, we use a 2.5- 
dimensional GRMHD code with Boyer-Lindquist coor- 
dinates {r,6,(j)). The method is based on a 3-1-1 for- 
malism of the general relativistic conservation laws of 
particle number and energy momentum, Maxwell equa- 
tions, and Ohm's law with no electrical resistance (ideal 
MHD condition) in a curved spacetime. In the simu- 
lations presented here we use minmod slope limiter re- 
construction, HLL approximate Riemann solver, flux- 
CT scheme and Noble's 2D method (see Mizuno et al. 
2006 and references therein). 

A geometrically thin Keplerian disk rotates around a 
black hole (non-rotating, a = 0.0 or rapidly co-rotating, 
a = 0.95, here a is black hole spin parameter), where 
the disk density is 100 times the coronal density. The 
thickness of the disk is H/r ^ 0.06. The background 
corona is free-falling, and the initial magnetic field is 
uniform and parallel to the rotational axis. Simula- 
tions are normalized by the speed of light, c, and the 
Schwarzschild radius, rs, with timescale, ts = rs/c. 
Values of the magnetic field strength and gas pres- 
sure depend on the normalized density, po- In these 
simulations the magnetic field strength, Bq, is set to 
0.05Vpoc2)- The 128 X 128 computational grid with 
logarithmic spacing in the radial direction spans the re- 
gion l.lrs < r < 20. Org (non-rotating black hole) and 
0.75rs < r < 20. Org (rapidly rotating black hole) and 
0.03 < 9 < 7r/2 where we assume axisymmetry with re- 
spect to the z-axis and mirror symmetry with respect to 
the equatorial plane. We employ a free boundary con- 
dition at the inner and outer boundaries in the radial 
direction. 

Figure 1 shows snapshots of the density (panels (a) 
and (b)), plasma beta (/? = Pgas/Pmag) distribution 
(panels (c) and (d)), and total velocity (panels (e) and 
(f)) for the non-rotating black hole, a = 0.0 (left pan- 
els); and the rapidly rotating black hole, a = 0.95 
(right panels); at each simulation's terminal time (non- 
rotating: t ~ 275Tg and rotating: t = 200Tg). At the 
marginally stable circular orbit (r = 3rs) the disk or- 
bits the black hole in about 40rs. The total velocity 



distribution of non-rotating and rapidly rotating black 
hole cases are shown in Figs, le and If. The jets in both 
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Fig. 1. — Snapshots of the non-rotating black hole (a, 
c, e) and the rapidly rotating black hole (6, d, f) at the 
applicable terminal simulation time. The color scales 
show the logarithm of density (upper panels), plasma 
beta (J5 — Pgas/Pmag; middle panels) and total velocity 
(lower panels). A negative velocity indicates infiow to- 
wards the black hole. The white lines indicate magnetic 
field lines (contour of the poloidal vector potential; up- 
per panels) and contours of the toroidal magnetic field 
strength (middle panels). Arrows depict the poloidal 
velocities normalized to light speed, as indicated above 
each panel by the arrow. 

cases have speeds greater than 0.4 c (mildly relativistic) 
that are comparable to the Alfven speeds. In the jets, 
toroidal velocity is the dominant velocity component. 
In the rapidly rotating black hole case, the velocity dis- 
tribution indicates a two-component jet with the inner 
jet not seen in the non-rotating black hole case. The 
inner jet is faster than the outer jet (over 0.5c). 
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3 RMHD Spine-Sheath Simulations 

In these simulations a "preexisting" jet is established 
across a computational domain of 6Rj x 6Rj x 60Rj with 
60 X 60 X 600 zones. The jet is in total pressure balance 
with a lower-density magnetized sheath with Pj/pe = 
2.0, where p is the mass density in the proper frame. 
The jet speed is Uj = 0.9165 c and 7^ = (1 — )~^/^ = 
2.5. The initial magnetic field is uniform and parallel 
to the jet flow. A precessional perturbation is applied 
at the inflow by imposing a transverse component of 
velocity with u± — O.Oluj. Here we show simulations 
with a precessional perturbation of angular frequency 
ojRj/uj = 0.93. In order to investigate the effect of 
an external wind, we have performed a no wind case 
{ue = 0) and a relativistic wind case (ue — 0.5 c). 
Simulations are halted after ~ 60 light crossing times 
of the jet radius (see Mizuno et al. 2007 for details). 

We have performed weakly magnetized simulations 
with sound speeds Ue ~ 0.57 c and Uj 0.51 c, and 
Alfven speeds vac ~ 0.07 c and VAj ~ 0.06 c. The sta- 
bilizing effect of a sheath wind is revealed in Figure 2. 
Here we see considerable reduction in transverse struc- 
k/ZD, no wind. a=0.9S. time=6().() 



RMHD, no wind, (x)=0.93, time=60.0 




Fig. 2. — 3D isovolume density image of the weakly 
magnetized case with no wind (top) and a c/2 wind 
(bottom). Magnetic field lines in white. 

ture and the jet spine reaches a larger distance before 
disruption in the presence of a wind. 

We have also performed strongly magnetized simula- 
tions with Alfven speeds VAe ~ 0.56 c and VAj ~ 0.45 c, 
and sound speeds Og 0.30 c and aj ~ 0.23 c. The 
stabilizing influence of a magnetic field and the stabi- 
lization of the jet spine in the presence of a magnetized 
sheath wind is shown in Figure 3. Here we see that the 
presence of the strong magnetic field has stabilized the 
jet spine even more than occured for the weakly mag- 
netized wind case and the initial helical perturbation 
is damped in the presence of the strongly magnetized 
sheath wind. 

More quantitatively we can analyse the growth or 
damping of the initial perturbation via ID cuts in the 




Fig. 3. — 3D isovolume density image of the strongly 
magnetized case with no wind (top) and a c/2 wind 
(bottom). Magnetic field lines in white. 
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Fig. 4. — Radial velocity (vx) along one dimensional 
cuts parallel to the jet axis and located at x/Rj — 0.2 
(solid line), 0.5 (dotted line) and 0.8 (dashed line) for 
the weakly magnetized (top) and strongly magnetized 
cases (bottom). 

radial velocity as shown in Figure 4. Measurable reduc- 
tion in transverse motion is seen for the weakly mag- 
netized wind case, significant reduction occurs for the 
strongly magnetized no wind case and stabilization oc- 
curs in the strongly magnetized wind case. 
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4 RMHD Spine-Sheath Stabihty 

Stability of a jet spine-sheath configuration can be an- 
alyzed by modeling the jet/spine as a cylinder of radius 
R embedded in an infinite sheath. A dispersion relation 
describing the growth or damping of the normal modes 
can be derived assuming uniform conditions within the 
spine, e.g., a uniform proper density, pj, axial magnetic 
field, Bj — Bj z, and velocity, — Uj^z, and assuming 
uniform conditions in the external sheath, e.g., a uni- 
form proper density, pe, axial magnetic field, B,, = B^^z, 
and velocity Ug — Ue,z (see Hardee 2007). 

Each normal mode consists of a single fundamental 
and multiple body wave solutions to the dispersion rela- 
tion. In the low frequency limit the helical fundamental 
mode has an analytical solution given by 
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In equation (2) ^Aj,e = (1 — v^j^^/c^) is the Alfven 
Lorentz factor. The jet is stable when 
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In the low frequency limit the real part of the first he- 
lical body mode has an analytical solution given by 
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Fig. 5. — Solutions for helical fundamental (red lines) 
and first body (green lines) modes for weakly magne- 
tized (flj^e ^ "^Aj^e) and strongly magnetized (wAj.e ~ 
2aj^e) jet simulations with no wind (ug = 0) and with a 
c/2 wind (u^ — 0.5). Solutions show the real, krRj, 
(dashed lines) and imaginary, kiRj, (dash-dot lines) 
parts of the wavenumber as a function of the angular 
frequency, ojRj/uj. Where the imaginary part of the 
wavenumber is shown in blue, the solution is damped. 
Immediately under the solutions for fundamental (S) 
and first body (Bl) modes is a panel that shows the 
wavelength, X/Rj, (dash-dot lines) and wave speed, 
v^/c, (dotted lines). The simulation precession fre- 
quency oj2 = 0.93 is indicated by the vertical solid line. 



In equation (4) Vms is a magnetosonic speed defined by 

1 /2 

Vms = ll\ + ''^a] where a is the sound speed, and 
VA is the Alfven wave speed. 

Equations (1 & 4) provide estimates for the helical 
fundamental and first body modes that can be followed 
by root finding techniques to higher frequencies. The 
results of numerical solution to the dispersion relation 
for the parameters appropriate to the numerical simu- 
lations shown in §3 are displayed in Figure 5. 

In the weakly magnetized cases fundamental (S) 
mode solutions consist of a growing (shown) and 
damped (not shown) solution pair (see eq. 1) and first 
body (Bl) mode solutions consist of a real and growing 
or damped solution pair. The presence of the external 
wind fiow leads to reduced growth of the S mode and 
weak damping of the Bl mode. 



In the strongly magnetized no wind case S mode so- 
lutions again consist of a growing and damped solution 
pair. However, we now find multiple growing solutions 
associated with the Bl mode at lower frequencies, and 
a modest damping rate accompanies the crossing of the 
multiple body mode solutions. At higher frequencies 
the Bl mode is similar to the weakly magnetized case. 

In the strongly magnetized wind case weak growth 
is associated with the slower, Ss^ moving shorter wave- 
length solution and weak damping is associated with 
the faster, 5/, moving longer wavelength solution. At 
frequencies, < lol, the growth rate is larger than the 
damping rate but at higher frequencies the damping 
rate is larger than growth rate for the S mode solution 
pair. In general the Bl mode is damped. 
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5 Conclusions 

Increased stability of the weakly-magnetized system 
with mildly relativistic sheath flow and stabilization of 
the strongly-magnetized system with mildly relativistic 
sheath flow is in agreement with theoretical results. In 
the fluid limit the present results conflrm earlier results 
obtained by Hardee & Hughes (2003), who found that 
the development of sheath flow around a relativistic jet 
spine explained the partial stabilization of the jets in 
their numerical simulations. 

The simulation results agree with theoretically pre- 
dicted wavelengths and wave speeds. On the other 
hand, growth rates and spatial growth lengths obtained 
from the linearized equations or from the present rela- 
tively low resolution simulations only provide guidelines 
to the rate at which perturbations grow or damp. 

A rapid decline in perturbation amplitudes in the 
sheath as a function of radius, governed by a Hankel 
function in the dispersion relation, suggests that the 
present results will apply to sheaths more than about 
three times the spine radius in thickness. 

Where flow and magentic flelds are parallel, current 
driven (CD) modes are stable (Isotomin & Pariev 1994, 
1996). However, we expect magnetic fields to have a 
significant toroidal component. Provided radial gradi- 
ents are not too large we expect the present results to 
remain approximately valid where Uj^e and Bje refer to 
poloidal velocity and magnetic field components. 

In the helically twisted magnetic and fiow field 
regime likely to be relevant to many astrophysical jets 
CD modes (Lyubarskii 1999) and/or KH modes could 
be unstable. While both CD and KH instability pro- 
duce helically twisted structure, the conditions for in- 
stability, the radial structure, the growth rate and the 
pattern motions are different. These differences may 
serve to identify the source of helical structure on rela- 
tivistic jet flows and allow determination of jet proper- 
ties near to the central engine. 
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